%% First run create_FRAP_structure.m to generate FRAP_struct for analysis
close all;
% run create_FRAP_structure.m; % Can comment this if already created structure
addpath(genpath("/Users/bknapp/Documents/FRAP_Data/analysis_scripts/"));

%% Plot each set of cap+cell fluor trace (as total fluorescence)
Ntraces = length(FRAP_struct);
nrows = floor(sqrt(Ntraces));
ncolumns = ceil(Ntraces/nrows);

figure;
for k=1:Ntraces
    timestamp = FRAP_struct(k).timestamp;
    time = timestamp - timestamp(1);
    cap_fluor = FRAP_struct(k).cap_fluor_tot;
    cell_fluor = FRAP_struct(k).cell_fluor_tot;
    
    subplot(nrows,ncolumns,k);
    plot(time,cell_fluor, 'k');
    hold on;
    plot(time,cap_fluor,'r');
    xlabel('Time (s)');
    ylabel('Fluor. (AU)');
end


%% Post-bleaching correction for photobleaching during recovery
frame_recovery = 3; % Frame at which recovery starts
fluor_background = 2.2; % Estimated mean pixel background intensity, may need to change this

figure;
for k=1:Ntraces
    % Obtain recovery data (raw)
    timestamp = FRAP_struct(k).timestamp;
    time_recovery = timestamp(frame_recovery:end) - timestamp(frame_recovery);
    cap_fluor_recovery = FRAP_struct(k).cap_fluor_tot(frame_recovery:end);
    cell_fluor_recovery = FRAP_struct(k).cell_fluor_tot(frame_recovery:end);
    
    % Generate background values
    cap_area = FRAP_struct(k).cap_area;
    cell_area = FRAP_struct(k).cell_area;
    
    cap_background = cap_area(frame_recovery:end)*fluor_background;
    cell_background = cell_area(frame_recovery:end)*fluor_background;
    
    % Correct for background intensity
    cell_fluor_recovery_bg = cell_fluor_recovery - cell_background;
    cap_fluor_recovery_bg = cap_fluor_recovery - cap_background;
    
    % Photobleaching effect can be estimated as fraction of total cell
    % fluor. that is photobleached during recovery (corrected for
    % background)
    FB = cell_fluor_recovery_bg./cell_fluor_recovery_bg(1);
    cell_fluor_recovery_bg_corrected = cell_fluor_recovery_bg./FB;
    cap_fluor_recovery_bg_corrected = cap_fluor_recovery_bg./FB;
    
    % Plot corrected recovery
    subplot(nrows,ncolumns,k);
    plot(time_recovery,cell_fluor_recovery_bg_corrected, 'k');
    hold on;
    plot(time_recovery,cap_fluor_recovery_bg_corrected,'r');
    xlabel('Time (s)');
    ylabel('Fluor. (AU)');
end

%% Plot FRAP recovery-only on single panel
frame_recovery = 3; % Frame at which recovery starts
fluor_background = 2.2; % Estimated mean pixel background intensity
mean_fluor_frames = 5; % Number of final frames to use for final recovered intensity

figure;
for k=1:Ntraces
    % Obtain recovery data (raw)
    timestamp = FRAP_struct(k).timestamp;
    time_recovery = timestamp(frame_recovery:end) - timestamp(frame_recovery);
    cap_fluor_recovery = FRAP_struct(k).cap_fluor_tot(frame_recovery:end);
    cell_fluor_recovery = FRAP_struct(k).cell_fluor_tot(frame_recovery:end);
    
    % Generate background values
    cap_area = FRAP_struct(k).cap_area;
    cell_area = FRAP_struct(k).cell_area;
    
    cap_background = cap_area(frame_recovery:end)*fluor_background;
    cell_background = cell_area(frame_recovery:end)*fluor_background;
    
    % Correct for background intensity
    cell_fluor_recovery_bg = cell_fluor_recovery - cell_background;
    cap_fluor_recovery_bg = cap_fluor_recovery - cap_background;
    
    % Photobleaching effect can be estimated as fraction of total cell
    % fluor. that is photobleached during recovery (corrected for
    % background)
    FB = cell_fluor_recovery_bg./cell_fluor_recovery_bg(1);
    cell_fluor_recovery_bg_corrected = cell_fluor_recovery_bg./FB;
    cap_fluor_recovery_bg_corrected = cap_fluor_recovery_bg./FB;
    
    % Calculate final corrected cap fluorescence value
    cap_final_value = mean(cap_fluor_recovery_bg_corrected(end-mean_fluor_frames:end));
    cap_initial_value = cap_fluor_recovery_bg_corrected(1);
    cap_fluor_diff = cap_final_value - cap_initial_value;
    cap_fluor_recovery_bg_corrected_norm = (cap_fluor_recovery_bg_corrected - cap_initial_value)/cap_fluor_diff;
    
    % Plot corrected recovery
    plot(time_recovery,cap_fluor_recovery_bg_corrected_norm, 'color', rand(1,3));
    hold on;
    xlabel('Time (s)');
    ylabel('Fluor. (AU)');
    
    % Perform fitting of exponential behavior
    fitfluor = 1 - cap_fluor_recovery_bg_corrected_norm;
    [f1,f2] = fit(time_recovery,fitfluor,'exp1');
    
    % Compute error in fit coefficients
    alpha = 0.95; % 95% confidence interval
    ci = confint(f1,alpha); % Obtain intervals of confidence (95%)
    df = length(time_recovery); % Degrees of freedom in T-test
    t = tinv((1+alpha)/2, df); % Use student's T test distribution
    se = (ci(2,:)-ci(1,:)) ./ (2*t); % Standard error of T-test
    
    % Now perform estimation using theoretical Gaussian spot for bleaching
    % (D = r^2/4*tauD); r = spot radius, tauD = characteristic recovery
    % time
    tauD = abs(1./f1.b); % Characteristic recovery time (min)
    tauD_err = abs(se(2)/f1.b^2);
    Deff = FRAP_struct(k).cap_diameter*FRAP_struct(k).cap_length/(4*tauD); % Effective diffusion coefficient (um^2/s)
    Deff_err = tauD_err/tauD^2;
    
    
    
    % Convert Deff to viscosity of membrane using brownian diffusion
    kB = 1.38E-23; % Boltzmann constant (J/K)
    T = 310; % Temperature (K)
    R = 2E-9; % Radius of probe molecule in meters (NOT SURE HOW GOOD THIS ESTIMATE IS)
    
    % Calculate effective viscosity based on Brownian diffusion model,
    % where D = kB*T/(6*pi*R*eta)
    eta_factor = kB*T/(6*pi*R);
    SI_conversion_Deff = (1E-12); % To convert Deff from um^2/s to m^2/s
    eta = eta_factor/(Deff*SI_conversion_Deff)*10; % Viscosity in Poise (P)
    eta_err = Deff_err/Deff^2;
    
    
    % Store data into FRAP_struct
    FRAP_struct(k).fluor_background = fluor_background;
    FRAP_struct(k).temperature = T;
    FRAP_struct(k).probe_radius = R;
    FRAP_struct(k).cap_final_fluor = cap_final_value;
    FRAP_struct(k).cap_fluor_recovery_corrected = cap_fluor_recovery_bg_corrected;
    FRAP_struct(k).time_recovery = time_recovery;
    FRAP_struct(k).Deff = Deff;
    FRAP_struct(k).Deff_err = Deff_err;
    FRAP_struct(k).eta = eta;
    FRAP_struct(k).eta_err = eta_err;
    
    
end
